Climate Change Forecasting Using Deep Learning

# import key libraries
import pandas as pd
import plotly.express as px
from copy import copy
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.preprocessing import MinMaxScaler
import tensorflow as tf
from tensorflow import keras
import plotly.offline as py
import plotly.express as px
import plotly.graph_objects as go

from jupyterthemes import jtplot
jtplot.style(theme = 'monokai', context = 'notebook', ticks = True, grid = False) 

# Data Source: https://www.kaggle.com/berkeleyearth/climate-change-earth-surface-temperature-data/notebooks?datasetId=29&sortBy=voteCount

# Read files
temperature_df = pd.read_csv('GlobalLandTemperaturesByCountry.csv')

temperature_df

temperature_df.shape

temperature_df.describe()

Exploratory data analysis

# Check the unique countries
temperature_df['Country'].unique()

# Check for missing values
temperature_df.isnull().sum()

# Check the dataframe information
temperature_df.info()

data cleaning

# Do groupby country to see the count of values available for each country
country_group_df = temperature_df.groupby(by = 'Country').count().reset_index('Country').rename(columns={'AverageTemperature':'AverageTemperatureCount','AverageTemperatureUncertainty' : 'AverageTemperatureUncertaintyCount'})

country_group_df

country_group_df['Country']

import plotly.express as px
fig = px.bar(country_group_df, x = 'Country', y = 'AverageTemperatureCount')
fig.show()

fig = px.bar(country_group_df, x = 'Country', y = 'AverageTemperatureUncertaintyCount')
fig.show()

# Plot histogram
fig = px.histogram(country_group_df, x = "AverageTemperatureCount")
fig.show()

fig = px.histogram(country_group_df, x = "AverageTemperatureUncertaintyCount")
fig.show()

country_group_df[(country_group_df['AverageTemperatureCount'] < 1500) | (country_group_df['AverageTemperatureUncertaintyCount'] < 1500)]

# Find countries with less than 1500 data info
countries_with_less_data = country_group_df[(country_group_df['AverageTemperatureCount'] < 1500) | (country_group_df['AverageTemperatureUncertaintyCount'] < 1500)].index.tolist()

countries_with_less_data

[8, 73, 80, 91, 98, 161, 204]

~temperature_df['Country'].isin(countries_with_less_data)

# Remove countries with less data info
temperature_df = temperature_df[~temperature_df['Country'].isin(countries_with_less_data)]

temperature_df.reset_index(inplace = True, drop = True)

temperature_df

# Fill missing values by doing rolling average on past 730 days
temperature_df['AverageTemperature'] = temperature_df['AverageTemperature'].fillna(temperature_df['AverageTemperature'].rolling(730, min_periods = 1).mean())

# Fill missing values by doing rolling average on past 730 days
temperature_df['AverageTemperatureUncertainty']= temperature_df['AverageTemperatureUncertainty'].fillna(temperature_df['AverageTemperatureUncertainty'].rolling(730, min_periods=1).mean())

temperature_df.isna().sum()

temperature_df['Country'].unique()

duplicates = []
for i in temperature_df['Country'].unique():
    if '(' in i:
        duplicates.append(i)
duplicates

# replace duplicates
temperature_df = temperature_df.replace(duplicates, ['Congo', 'Denmark','Falkland Islands','France','Netherlands','United Kingdom'])

temperature_df['Country'].unique()

data visualization

countries = temperature_df['Country'].unique().tolist()
countries

# mean temperature for each country
mean_temperature = []
for i in countries:
    mean_temperature.append(temperature_df[temperature_df['Country'] == i]['AverageTemperature'].mean())

# Plot mean teamperature of countries
data = [ dict(type = 'choropleth', # type of map
              locations = countries, # location names
              z = mean_temperature, # temperature of countries
              locationmode = 'country names')
       ]

layout = dict(title = 'Average Global Land Temperatures',
              geo = dict(showframe = False,
                         showocean = True, # to show the ocean
                         oceancolor = 'aqua',
                         projection = dict(type = 'orthographic'))) # to get the globe view),

fig = dict(data = data, layout = layout)
py.iplot(fig, validate = False, filename = 'worldmap')

# year of recorded data for visualization 
temperature_df['year'] = temperature_df['dt'].apply(lambda x: x.split('-')[0])
temperature_df

# creating the animation to see the global temperature change
fig = px.choropleth(temperature_df, locations = 'Country',
                    locationmode = 'country names', # locations 
                    color = 'AverageTemperature', # column representing the temperature
                    hover_name = "Country", # column to add to hover information
                    animation_frame = 'year', # timeframe for animation
                    color_continuous_scale = px.colors.sequential.deep_r)
# py.plot(fig)
fig.show()

# To get  global average tempeature over years
df_global = temperature_df.groupby('year').mean().reset_index()

df_global['year'] = df_global['year'].apply(lambda x: int(x))
df_global = df_global[df_global['year'] > 1850]

# uncertainity upper bound 
trace1 = go.Scatter(
    x = df_global['year'], 
    y = np.array(df_global['AverageTemperature']) + np.array(df_global['AverageTemperatureUncertainty']), # Adding uncertinity
    name = 'Uncertainty top',
    line = dict(color = 'green'))

# uncertainity lower bound
trace2 = go.Scatter(
    x = df_global['year'] , 
    y = np.array(df_global['AverageTemperature']) - np.array(df_global['AverageTemperatureUncertainty']), # Subtracting uncertinity
    fill = 'tonexty',
    name = 'Uncertainty bottom',
    line = dict(color = 'green'))

# recorded temperature
trace3 = go.Scatter(
    x = df_global['year'] , 
    y = df_global['AverageTemperature'],
    name = 'Average Temperature',
    line = dict(color='red'))
data = [trace1, trace2, trace3]

layout = go.Layout(
    xaxis = dict(title = 'year'),
    yaxis = dict(title = 'Average Temperature, °C'),
    title = 'Average Land Temperatures Globally',
    showlegend = False)

fig = go.Figure(data = data, layout = layout)
py.iplot(fig)

# modeling for US data
US_df = temperature_df[temperature_df['Country'] == 'United States'].reset_index(drop = True)
US_df

fig = px.line(title = 'US Temperature Data')
US_df_updated = US_df[US_df['year'] > '1813']
fig.add_scatter(x = US_df_updated['dt'], y = US_df_updated['AverageTemperature'], name = 'US Temperature')
fig.show()

data preparation for model training

# Get the month of recording, to use as a feature
temperature_df['Month'] = temperature_df['dt'].apply(lambda x: int(x.split('-')[1]))

# To get the global average tempeature over years
df_global_monthly = temperature_df.groupby(['dt']).mean().reset_index()

df_global_monthly

# creating data for training the time series model
def prepare_data(df, feature_range):
    # Get the columns
    columns = df.columns
    # For the given range, create lagged input feature for the given columns
    for i in range(1, (feature_range + 1)):
        for j in columns[1:]:
            name = j + '_t-' + str(i)
            df[name] = df[j].shift((i))
    # Create the target by using next value as the target
    df['Target'] = df['AverageTemperature'].shift(-1)
    return df

df_global_monthly = prepare_data(df_global_monthly, 3)

df_global_monthly

df_global_monthly = df_global_monthly.dropna().reset_index(drop = True)

df_global_monthly

train = df_global_monthly[:int(0.9 * len(df_global_monthly))].drop(columns = 'dt').values

test = df_global_monthly[int(0.9 * len(df_global_monthly)):].drop(columns = 'dt').values

scaler = MinMaxScaler(feature_range = (0, 1))
train  = scaler.fit_transform(train)
test   = scaler.transform(test)

# split data into input features, targets
train_x, train_y = train[:,:-1], train[:,-1]
test_x, test_y = test[:,:-1], test[:,-1]

# reshape input to 3D [samples, timesteps, features]
train_x = train_x.reshape((train_x.shape[0], 1, train_x.shape[1]))
test_x = test_x.reshape((test_x.shape[0], 1, test_x.shape[1]))
print(train_x.shape, train_y.shape, test_x.shape, test_y.shape)

(2911, 1, 12) (2911,) (324, 1, 12) (324,)

LSTM Model for predicting global temperature trend

def create_model(train_x):
    # create model
    inputs = keras.layers.Input(shape = (train_x.shape[1], train_x.shape[2]))
    x = keras.layers.LSTM(50,return_sequences =  True)(inputs)
    x = keras.layers.Dropout(0.3)(x)
    x = keras.layers.LSTM(50, return_sequences = True)(x)
    x = keras.layers.Dropout(0.3)(x)
    x = keras.layers.LSTM(50)(x)
    outputs = keras.layers.Dense(1, activation = 'linear')(x)

    model = keras.Model(inputs = inputs, outputs = outputs)
    model.compile(optimizer = 'adam', loss = "mse")
    return model

model = create_model(train_x)
model.summary()

Model: "model"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_1 (InputLayer)         [(None, 1, 12)]           0         
_________________________________________________________________
lstm (LSTM)                  (None, 1, 50)             12600     
_________________________________________________________________
dropout (Dropout)            (None, 1, 50)             0         
_________________________________________________________________
lstm_1 (LSTM)                (None, 1, 50)             20200     
_________________________________________________________________
dropout_1 (Dropout)          (None, 1, 50)             0         
_________________________________________________________________
lstm_2 (LSTM)                (None, 50)                20200     
_________________________________________________________________
dense (Dense)                (None, 1)                 51        
=================================================================
Total params: 53,051
Trainable params: 53,051
Non-trainable params: 0
_________________________________________________________________

# fit the network
history = model.fit(train_x, train_y, epochs = 50, batch_size = 72, validation_data = (test_x, test_y), shuffle = False)

Train on 2911 samples, validate on 324 samples
Epoch 1/50
2911/2911 [==============================] - 10s 3ms/sample - loss: 0.3119 - val_loss: 0.0426
Epoch 2/50
2911/2911 [==============================] - 0s 165us/sample - loss: 0.0847 - val_loss: 0.0289
Epoch 3/50
2911/2911 [==============================] - 0s 168us/sample - loss: 0.0579 - val_loss: 0.0185
Epoch 4/50
2911/2911 [==============================] - 0s 164us/sample - loss: 0.0418 - val_loss: 0.0107
Epoch 5/50
2911/2911 [==============================] - 0s 137us/sample - loss: 0.0299 - val_loss: 0.0038
Epoch 6/50
2911/2911 [==============================] - 1s 184us/sample - loss: 0.0190 - val_loss: 0.0020
Epoch 7/50
2911/2911 [==============================] - 0s 170us/sample - loss: 0.0155 - val_loss: 0.0017
Epoch 8/50
2911/2911 [==============================] - 0s 169us/sample - loss: 0.0133 - val_loss: 0.0016
Epoch 9/50
2911/2911 [==============================] - 1s 188us/sample - loss: 0.0132 - val_loss: 0.0016
Epoch 10/50
2911/2911 [==============================] - 1s 179us/sample - loss: 0.0124 - val_loss: 0.0017
Epoch 11/50
2911/2911 [==============================] - 0s 139us/sample - loss: 0.0117 - val_loss: 0.0014
Epoch 12/50
2911/2911 [==============================] - 0s 150us/sample - loss: 0.0107 - val_loss: 0.0014
Epoch 13/50
2911/2911 [==============================] - 0s 154us/sample - loss: 0.0107 - val_loss: 0.0014
Epoch 14/50
2911/2911 [==============================] - 1s 237us/sample - loss: 0.0100 - val_loss: 0.0014
Epoch 15/50
2911/2911 [==============================] - 0s 149us/sample - loss: 0.0100 - val_loss: 0.0013
Epoch 16/50
2911/2911 [==============================] - 0s 137us/sample - loss: 0.0092 - val_loss: 0.0012
Epoch 17/50
2911/2911 [==============================] - 0s 150us/sample - loss: 0.0093 - val_loss: 0.0012
Epoch 18/50
2911/2911 [==============================] - 1s 207us/sample - loss: 0.0088 - val_loss: 0.0011
Epoch 19/50
2911/2911 [==============================] - 1s 184us/sample - loss: 0.0085 - val_loss: 9.9897e-04
Epoch 20/50
2911/2911 [==============================] - 1s 200us/sample - loss: 0.0079 - val_loss: 0.0010
Epoch 21/50
2911/2911 [==============================] - 1s 196us/sample - loss: 0.0082 - val_loss: 9.4213e-04
Epoch 22/50
2911/2911 [==============================] - 0s 171us/sample - loss: 0.0081 - val_loss: 0.0010
Epoch 23/50
2911/2911 [==============================] - 1s 173us/sample - loss: 0.0078 - val_loss: 0.0010
Epoch 24/50
2911/2911 [==============================] - 0s 162us/sample - loss: 0.0079 - val_loss: 0.0011
Epoch 25/50
2911/2911 [==============================] - 0s 171us/sample - loss: 0.0074 - val_loss: 0.0011
Epoch 26/50
2911/2911 [==============================] - 1s 179us/sample - loss: 0.0076 - val_loss: 9.5496e-04
Epoch 27/50
2911/2911 [==============================] - 1s 195us/sample - loss: 0.0074 - val_loss: 0.0010
Epoch 28/50
2911/2911 [==============================] - 0s 157us/sample - loss: 0.0073 - val_loss: 0.0011
Epoch 29/50
2911/2911 [==============================] - 1s 202us/sample - loss: 0.0072 - val_loss: 0.0011
Epoch 30/50
2911/2911 [==============================] - 1s 265us/sample - loss: 0.0070 - val_loss: 0.0011
Epoch 31/50
2911/2911 [==============================] - 1s 257us/sample - loss: 0.0068 - val_loss: 0.0010
Epoch 32/50
2911/2911 [==============================] - 1s 222us/sample - loss: 0.0071 - val_loss: 0.0011
Epoch 33/50
2911/2911 [==============================] - 1s 211us/sample - loss: 0.0070 - val_loss: 0.0011
Epoch 34/50
2911/2911 [==============================] - 1s 208us/sample - loss: 0.0069 - val_loss: 0.0011
Epoch 35/50
2911/2911 [==============================] - 1s 194us/sample - loss: 0.0069 - val_loss: 0.0011
Epoch 36/50
1368/2911 [=============>................] - ETA: 0s - loss: 0.0101

def plot_history(history):
    # plot history
    plt.plot(history.history['loss'], label='train')
    plt.plot(history.history['val_loss'], label='test')
    plt.grid()
    plt.legend()
    plt.show()

plot_history(history)

# model performance 
def prediction(model,test_x,train_x, df):
    # Predict using the model
    predict =  model.predict(test_x)

    # Reshape test_x and train_x for visualization  and inverse-scaling purpose
    test_x = test_x.reshape((test_x.shape[0], test_x.shape[2]))
    train_x = train_x.reshape((train_x.shape[0], train_x.shape[2]))

    # Concatenate test_x with predicted value
    predict_ = np.concatenate((test_x, predict),axis = 1)

    predict_ = scaler.inverse_transform(predict_)
    original_ = scaler.inverse_transform(test)

    # Create dataframe to store the predicted and original values
    pred = pd.DataFrame()
    pred['dt'] = df['dt'][-test_x.shape[0]:]
    pred['Original'] = original_[:,-1]
    pred['Predicted'] = predict_[:,-1]

    pred['Error'] = pred['Original'] - pred['Predicted']
    
    # Create dataframe for visualization
    df = df[['dt','AverageTemperature']][:-test_x.shape[0]]
    df.columns = ['dt','Original']
    original = df.append(pred[['dt','Original']])
    df.columns = ['dt','Predicted']
    predicted = df.append(pred[['dt','Predicted']])
    original = original.merge(predicted, left_on = 'dt',right_on = 'dt')
    return pred, original

pred, original = prediction(model, test_x, train_x, df_global_monthly )

def plot_error(df):

    # Plotting the Current and Predicted values
    fig = px.line(title = 'Prediction vs. Actual')
    fig.add_scatter(x = df['dt'], y = df['Original'], name = 'Original', opacity = 0.7)
    fig.add_scatter(x = df['dt'], y = df['Predicted'], name = 'Predicted', opacity = 0.5)
    fig.show()

    fig = px.line(title = 'Error')
    fig = fig.add_scatter(x = df['dt'], y = df['Error'])
    fig.show()

def plot(df):
    # Plotting the Current and Predicted values
    fig = px.line(title = 'Prediction vs. Actual')
    fig.add_scatter(x = df['dt'], y = df['Original'], name = 'Original', opacity = 0.7)
    fig.add_scatter(x = df['dt'], y = df['Predicted'], name = 'Predicted', opacity = 0.5)
    fig.show()

plot(original)

plot_error(pred)

data prep for US data

US_df = temperature_df[temperature_df['Country'] == 'United States'].reset_index(drop = True)
US_df

US_df = US_df.drop(['Country', 'year'], axis = 1)

# training data
US_df = prepare_data(US_df, 3)

US_df

US_df = US_df.dropna().reset_index(drop = True)

train = US_df[:int(0.9 * len(US_df))].drop(columns = 'dt').values

test = US_df[int(0.9 * len(US_df)):].drop(columns = 'dt').values

train.shape

scaler = MinMaxScaler(feature_range = (0, 1))
train  = scaler.fit_transform(train)
test   = scaler.transform(test)

# split data into input features and targets
train_x, train_y = train[:,:-1], train[:,-1]
test_x, test_y = test[:,:-1], test[:,-1]

# reshape input to 3D 
train_x = train_x.reshape((train_x.shape[0], 1, train_x.shape[1]))
test_x = test_x.reshape((test_x.shape[0], 1, test_x.shape[1]))
print(train_x.shape, train_y.shape, test_x.shape, test_y.shape)

train_x

LSTM Model for predicting US temperature trend

model = create_model(train_x)
model.summary()

# fit the network
history = model.fit(train_x, train_y, epochs = 50, batch_size = 72, validation_data = (test_x, test_y), shuffle = False)

plot_history(history)

# us model performance 
pred, original = prediction(model, test_x, train_x, US_df )

plot(original)

plot_error(pred)